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A calculation of the shear viscosity in SU(3) gluodynamics 
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We perform a lattice Monte-Carlo calculation of the two-point functions of the energy-momentum 
tensor at finite temperature in the SU(3) gauge theory. Unprecedented precision is obtained thanks 
to a multi-level algorithm. The lattice operators are renormalized non-perturbatively and the clas- 
sical discretization errors affecting the correlators are corrected for. A robust upper bound for the 
shear viscosity to entropy density ratio is derived, rj/s < 1.0, and our best estimate is rj/s = 0.1 34(33) 
at T = 1.65T C under the assumption of smoothness of the spectral function in the low- frequency 
region. 

PACS numbers: 12.38.Gc, 12.38.Mh, 25.75.-q 



Introduction. — Models treating the system produced in 
heavy ion collisions at RHIC as an ideal fluid have had 
significant success in describing the observed flow phe- 
nomena [l|, Q . Subsequently the leading corrections due 
to a finite shear viscosity were computed [3j, in parti- 
cular the flattening of the elliptic flow coefficient V2(pt) 
above lGeV. It is therefore important to compute the 
QCD shear and bulk viscosities from first principles to 
establish this description more firmly. Small transport 
coefficients are a signature of strong interactions, which 
lead to efficient transmission of momentum in the system. 
Strong interactions in turn require non-perturbative com- 
putational techniques. Several attempts have been made 
to compute these observables on the lattice in the SU(3) 
gauge theory 0, @ • The underlying basis of these calcu- 
lations are the Kubo formulas, which relate each trans- 
port coefficient to a spectral function p(u) at vanishing 
frequency. Even on current computers, these calcula- 
tions are highly non-trivial, due to the fall-off of the re- 
levant correlators in Euclidean time (as Xq 5 at short dis- 
tances), implying a poor signal-to-noise ratio in a stan- 
dard Monte-Carlo calculation. The second difficulty is 
to solve the ill-posed inverse problem for p(ui) given the 
Euclidean correlator at a finite set of points. Mathemati- 
cally speaking, the uncertainty on a transport coefficient 
X is infinite for any finite statistical accuracy, because 
adding eujS(u>) to p(u>) merely corresponds to adding a 
constant to the Euclidean correlator of order e, while ren- 
dering x infinite. Therefore smoothness assumptions on 
p{u>) have to be made, which are reasonable far from the 
one-particle energy eigenstates, and can be proved in the 
hard-thermal- loop framework [6| . 

In this Letter we present a new calculation which dra- 
matically improves on the statistical accuracy of the Eu- 
clidean correlator relevant to the shear viscosity through 
the use of a two-level algorithm [16| . This allows us to 
derive a robust upper bound on the viscosity and a use- 



ful estimate of the ratio rj/s, which has acquired a special 
significance since its value 1 /4ir in a class of strongly cou- 
pled supersymmetric gauge theories Q was conjectured 
to be an absolute lower bound for all substances @ . 
Methodology. — In the continuum, the energy- 
momentum tensor T^{x) = F^ a F^ a - \5^F^F^ a , be- 
ing a set of Noether currents associated with translations 
in space and time, does not renormalize. With Lq = 1/T 
the inverse temperature, we consider the Euclidean two- 
point function (0 < xq < Lq) 



C(x ) = Ll 



(T 12 (0)T 12 (z ,x)). 



(1) 



The tree-level expression is C t L (xo) = {f{r) - ^) , 
with t = 1 — ^a, d A = 8 the number of gluons and 

f(z) = ds s 4 cosh 2 (zs)/ sinh 2 s. The correlator C(xq) 
is thus dimensionless and, in a conformal field theory, 
would be a function of Txq only. 
The spectral function is defined by 



C(x ) = Ll / p(u) 



t coshw(i-Lo — ^o) 
sinh^a 



du>. 



The shear viscosity is given by 0, 0j| 



T)(T) = 7T -f- 

du> 



(2) 



(3) 



Important properties of p are its positivity, p(u))/u) > 
and parity, p(—w) = —p(u>). The spectral function that 
reproduces C t,L (iq) is 



A tx = 



Ax ^ 
tanh juiL 

1 d A 
10 (4tt) 2 ' 



+ BL 4 lj6(u), 



15 



d A - 



(4) 
(5) 
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While the w 4 term is expected to survive in the inter- 
acting theory with only logarithmic corrections, the 5- 
function at the origin corresponds to the fact that gluons 



2 



are asymptotic states in the free theory and implies an 
infinite viscosity. 

On the lattice, translations only form a discrete group, 
so that a finite renormalization is necessary, T liv {go) = 

Z(go)T^ Te . We employ the Wilson action [9], S s = 
^E I , /I ^ Tr i 1 - P^{x)}, on an L ■ L 3 hypertoroidal 
lattice, and the following discretized expression of the 
Euclidean energy: 



T, 



(bare) 
00 



(x) 



« 4 5o 2 



[^ReTrP fei (x) 

k<l 



^ReTrP 0fe (a 



One of the lattice sum rules [l(| can be interpreted as a 
non-perturbative renormalization condition for this par- 
ticular discretization, from which we read off Z(go) = 
1 — \g^(c a — c r ). The definition of the anisotropy coeffi- 
cients c CTjT can be found in [12|, where they are computed 
non-perturbatively. With a precision of about 1%, a Pade 
fit constrained by the one-loop result [ll| yields 



Z(go) = 



1 - 1.0225^ + 0.1305^ 
1 - 0.8557c/ 2 : 



(6/g 2 > 5.7). (6) 



Numerical results. — We report results obtained on a 
(3 = 6/g% = 6.2, 8 • 20 3 lattice and on a (3 = 6.408, 8 • 28 3 
lattice. The first is thus at a temperature of 1.24T C , the 
second at T = 1.65T C . We use the results for aT c obtained 
in [3| and the non-perturbative lattice /3-function of [Hj] 
to determine this. We employ the two-level algorithm 
described in [Til ]. The computing time invested into the 
1.65T C simulation is about 860 PC days. Following we 
discretize i((T u - T 22 )(T U - T 22 )) instead of (T 12 T 12 ) 
(the two are equal in the continuum) to write C(xq) = 



%{O v (0)O n (x )) + 0(a 2 ), where 

O v (x Q ) = ±a 3 ^{T n -T 22 }( 5o ,: 
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FIG. 1: The correlators that contribute to C(xo) = \{Cbb + 
Gee + ICeb). Filled symbols correspond to T = 1.65T C , 
open symbols to 1.24T C . Error bars are smaller than the data 
symbols. 
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The three electric-electric, magnetic-magnetic and 
electric-magnetic contributions to C(xq) are computed 
separately and shown on Fig. Q] We apply the follow- 
ing technique to remove the tree-level discretization er- 
rors [ljj separately to Cbb, Cee and Ceb- Firstly, xq 
is defined such that C^ nt (xo) = C^(x ). The improved 
correlator is defined at a discrete set of points through 
C(xq) = C(xq), and then augmented to a continuous 
function via C{x^) = a + pC^ nt (x^), i = 1,2, where 

Xq and Xq 2 '' correspond to two adjacent measurements. 

The resulting improved correlator, normalized by the 
continuum tree-level result, is shown on Fig. [2l One ob- 
serves that the deviations from the tree-level result are 
surprisingly small, while deviations from conformality are 
visible. The latter is not unexpected at these tempera- 
tures, where p/T 4 is still strongly rising [l8|]. Finite- 
volume effects on the T — 1.65T C lattice are smaller than 



FIG. 2: The tree-level improved correlator C(xo) normalized 
to the tree-level continuum infinite-volume prediction. The 
four points in each sequence are strongly correlated, but their 
covariance matrix is non-singular. 



one part in 10 3 at tree-level. Non-perturbatively, at the 
same temperature with resolution L$/a — 6, increasing 
L/a from 20 to 30 reduces C(L /2) by a factor 0.922(73). 
While not statistically compelling as it stands, the effect 
deserves further investigation. 

The entropy density is obtained from the relation 
s = (e + p)/T and the standard method to compute 
e + p (12], Eq.1.14). We find s/T 3 = 4.72(3)(5) and 
5.70(2)(6) respectively at T/T c = 1.24 and 1.65 (the first 
error is statistical and the second is the uncertainty on 
Z(go)). The Stefan-Boltzmann value is 32-7r 2 /45 in the 
continuum and 1.0867 times that value [3] at Lo/a = 8. 
Unsatisfactory attempts to extract the viscosity. — In 
order to compare with previous studies 0,@,wefit C(x ) 
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with a Breit-Wigner ansatz 



F 



F 



l + b 2 (LU~LU ) 2 l + b 2 (uj + LU a ) 2 



(7) 



although it clearly ignores asymptotic freedom, which im- 
plies that p(uj) ~ w 4 at w > T [|. The result of a cor- 
related fit at T = 1.65T C using the points at Txq = 0.5, 
0.35 and 0.275 is a 3 F = 0.78(4), (b/a) 2 = 240(30) and 
au} = 2.36(4), and hence r) / s\t=i.65T <1 _= 0.33(3). A 
comparison of this to the results of Ref. Q illustrates the 
progress made in statistical accuracy. 

An ansatz motivated by the hard-thermal-loop frame- 
work is m 



p(w)/u 



r)/ir 



l + b 2 u 



+ 9(w — LO\) 



Alu 3 



tanhw/4T' 



(8) 



It is capable of reproducing the tree-level prediction, 
Eq. 01 and it allows for a thermal broadening of the delta 
function at the origin. Fitting the T — 1.65T C points 
shown on Fig. ^ the x 2 is minimized for b — (effec- 
tively eliminating a free parameter), A/An, — 0.996(8), 
Wi/T = 7.5(2) and rj/s = 0.25(3), with x^n = 4 -°- Thus 
while the ansatz is hardly compatible with the data, it 
shows that the data tightly constrains the coefficient A 
to assume its tree-level value. 

A bound on the viscosity. — The positivity property of 
p(u)) allows us to derive an upper bound on the viscosity, 
based on the following assumptions: 

1. the contribution to the correlator from u> > A is 
correctly predicted by the tree-level formula 

2. the width of any potential peak in the region to < T 
is no less than O(T). 

The standard QCD sum rule practice is to use perturba- 
tion theory from the energy lying midway between the 
lightest state and the first excitation. With this in mind 
we choose A = max(±[M 2 + M 2 »] » 2.6GeV, 5T), where 
M 2 (») are the masses of the two lightest tensor glueballs. 
Perturbation theory predicts a Breit-Wigner centered at 
the origin of width T = 2 7 @, where 7 ~ a s NT is the 
gluon damping rate. To derive the upper bound we con- 
servatively assume that for lj < \/2T, p(uS)/oj is a Breit- 
Wigner of width r = T centered at the origin. From 



J PBW{UJ) + J A 



Pt.l.M 



sinhwLo/2 



we 



obtain (with 90% statistical confidence level) 



v/s< 



0.96 
1.08 



(T 
(T 



1.65T C 
1.24T r 



(9) 



The spectral function. — As illustrated above, it is rather 
difficult to find a functional form for p(uj) that is both 
physically motivated and fits the data. In a more model- 
independent approach, p(uj) is expanded in an orthogo- 
nal set of functions, which grows as the lattice resolution 
on the correlator increases, and becomes complete in the 




FIG. 3: The result for p(u>). The meaning of the error bands 
and the curves is described in the text. The area under them 
equals C{L /2) = 8.05(31) and 9.35(42) for 1.24T C and 1.65T C 
respectively. 



limit of Lq/o — > 00. We proceed to determine the func- 
tion p{uj) = p(w) / tanh(icjL ) by making the ansatz 



p(uj) = m(uj) [1 4- a(ui)], 



(10) 



where m(uj) > has the high-frequency behavior of 
Eq.01 and correspondingly define K(xq, uj) = coshci;(a;o — 
\Lq)/ cosh \uLq. Suppose that m(uj) already is a 
smooth approximate solution to p(u>); inserting (fTO]) into 
(Eq. [2|), one requires that a(uj) = ^^c^a^(w), with 
{ai} a basis of functions which is as sensitive as pos- 
sible to the discrepancy between the lattice correla- 
tor and the correlator generated by m(uj). These are 
the eigcnfunctions of largest eigenvalue of the symmet- 
ric kernel G(w,a/) = Jo° ^M(x ,uj)M(x ,uj'), where 
M(xq,u>) = K(xQ,uj)m(Lu). These functions satisfy 
J Q dajui(u>)ui> (lu) = Su' and have an increasing number 
of nodes as their eigenvalue decreases. Thus the more 
data points available, the larger the basis and the finer 
details of the spectral function one is able to determine. 

To determine the spectral function from N points of 
the correlator, we proceed by first discretizing the u> vari- 
able into an A^-vector. The final spectral function is 
given by the last member p( N ~> of a sequence whose first 
member is p^ — m and whose general member p( n ' re- 
produces n points (or linear combinations) of the lattice 
correlator. For n > 1, p^ = ^""^[l + £™=i cf^af ] ] 
and the functions a^{uj) are found by the SVD decom- 
position [H[ of the JV.xn matrix M^, where M^f = 

(luj). The 'model' m(uj) is thus up- 
dated and agrees with p(uj) at the end of the procedure. 
We first performed this procedure on coarser lattices 
with Lq/o = 6 at the same temperatures, starting from 
m(oj) = A t .i.ti' 4 /(tanh(-|wLo) tanh(iwLo) tanh 2 (cwLo)) 
with j < c < i , and then recycled the output as seed for 
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the Lo/a — 8 lattices. On the latter we used the N = 4 
points shown on Fig. [21 

The next question to address is the uncertainty on 
p(u>). It is important to realize that even in the absence 
of statistical errors, a systematic uncertainty subsists due 
to the finite number of basis functions we can afford to 
describe p{to) with. A reasonable measure of this uncer- 
tainty is by how much varies if one doubles the re- 
solution on C(xq). This can be estimated by 'generating' 
new points by using the computed p^ N ^(uj). On the other 
hand we perform a two-point interpolation in x -space 
(we chose the form (a + P(xo — ^L ) 2 )/ sin 5 (7nro/Lo))> 
and take the difference between these and the generated 
ones as their systematic uncertainty. In practice this dif- 
ference is added in quadrature with the statistical uncer- 
tainty. Next we repeat the procedure to find p described 
above with N — > 2N: if we use as seed p^ N \ then by 
construction it is left invariant by the iterative proce- 
dure, but the derivatives of p( 2JV ) with respect to the 
2N points of the correlator can be evaluated. The er- 
ror on p{oj) is then obtained from a formula of the type 

(Sp) 2 = J2l=i(~&c~) 2 (3Ci) 2 which however keeps track of 
correlations in x and Monte-Carlo time. This is the er- 
ror band shown on Fig. [3] and the corresponding shear 
viscosity values are 

, _ f 0.134(33) (T= 1.65T C ) , * 

v/s \ 0.102(56) (T=1.24T C ). [ ' 

It is also interesting to check for the stability of the so- 
lution under the use of a larger basis of functions. If in- 
stead of starting from p( N >(u>) we restart from p^ (the 
output of the L /a = 6 lattice) and fit the 2N (depen- 
dent) points using 27V basis functions {ai }, we obtain the 



curves drawn on Fig. [3l As one would hope, the oscilla- 
tions of p( 2Ar ) (oj) are covered by the error band. 

Conclusion. — Using state-of-the-art lattice techniques, 
we have computed the correlation functions of the energy- 
momentum tensor to high accuracy in the SU(3) pure 
gauge theory. We have calculated the leading high- 
temperature cutoff effects and removed them from the 
correlator relevant to the shear viscosity, and we nor- 
malized it non-perturbatively, exploiting existing results. 
We obtained the entropy density with an accuracy of 1%. 
The most robust result obtained on the shear viscosity 
is the upper bound Eq. ([9]), which comes from lumping 
the area under the curve on Fig. [3] in the interval [0, 6T] 
into a peak of width T = T centered at the origin. Sec- 
ondly, our best estimate of the shear viscosity is given by 
Eq. (flT|) , using a new method of extraction of the spectral 
function. The errors contain an estimate of the systema- 
tic uncertainty associated with the limited resolution in 
Euclidean time. We are extending the calculation to finer 
lattice spacings and larger volumes to further consolidate 
our findings. 

The values dTl"]) are intriguingly close to saturating the 
KSS bound 8] r//s > 1/Att. We note that in perturba- 
tion theory the ratio rj/s does not depend strongly on the 
number of quark flavors [20I ] . Our results thus corrobo- 
rate the picture of a near-perfect fluid that has emerged 
from the RIIIC experiments, with the magnitude of the 
anisotropic flow incompatible with rj/s > 0.2 [3J]. 
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